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** Abstract : The Lame's differential equation is solved by the finite-difference method. 
** Subjects: Viscoelasticity , Hereditary Solid Mechanics , The Differential equation . 



NOTE: This worksheet demonstrates the use of Maple for calculating the solution of Lame's differential 

equation . 

The authors expect that this worksheet will only be used for teaching and educational purposes .. 
Copyright 

Co.H Tran - Phong T Ngo - FINITE DIFFERENCE METHOD AND 

THE LAME'S EQUATION IN HEREDITARY SOLID MECHANICS . Use MliplC 9.5 
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A. THE DISPLACEMENT DIFFERENTIAL 
EQUATION : 



The Lame's equation of the plane-deformation problem in the cylinder made of orthotropic viscoelastic 
composite material does not have constant modules . 

The modules Et , Er and Ert will be replaced with the functions Er(t) , Et(t) and Ert(t) respectively . 
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a*. The plane-deformation problem : Bai toan bien dang phang cua ong tru true huong composite dan 
nhot : 

We examine an orthotropic viscoelastic composite material cylinder which has the horizontal section within 

limit of 2 circles :r = a,r = b(a<b). Choosing the cylindrical coordinates r , ^ , z ( the 
axial z is along with the cylinder ) . The components of stress and deformation are functions of r , t 
respectively . 

Xet ong tru co tiet dien ngang gioi han boi 2 duong tron dong tam co ban kinh r = a,r=b ( a < b ) , 
ong tru duoc lam bang vat lieu co tfnh true huong . Chon he toa do tru r , p , z ( true z huong doc 

theo ong tru . Cac thanh phan bien dang va ung suat tuong ung la ^" ' ^ e ' ^ r ' ^ e la cac ham theo r 
, t . 



The two components of deformation-tensor : ( 2 thanh phan cua tensor bien dang la : ) 

, . du(r,t) u(r,/) 
£ T (r,0 = * ; £,(rj) = 

or r 
and the differential equation of equilibrium is : ( phuong trinh vi phan can bang ) 

dajr.t) 

® r when t = 0 , boundary conditions : ( khi t = 0 , cac dieu kien 

bien : ) ^,(a,0) =-P ; a,(bfi) = -Q 

b* . The displacement - differential equation : Phuong trinh vi phan chuyen vi : 

The differential equation of the cylinder displacement in the case of viscoelastic plane-deformation : ( 
Phuong trinh chuyen vi ong trong truong hop bien dang phang dan nhot ) 

r = — + t ^— u(r. f) = 0 

dr dr E T 
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( To : const ) 



B. FINITE DIFFERENCE METHOD : 

The boundary conditions of the problem are given at two edges ( Dieu kien bien cua bai toan duoc cho o 2 
canh ) : r = a and r = b . (a=l,b = 2) 

Now we choose the number of mesh points ( Ta chon so diem luoi ) N = 20 . The interval over which we 
approximate this equation is ( Doan xap xi cua phuong trinh la ) [a, b] . And the step size for this interval is 

1 

( Va kich thuoc buoc nhay cho doan nay la ) 



The difference operators are ( Cac toan tu sai phan la ) Uj and Ujj , And we have two boundary 

:= t = 1 := lijo t = -1 

conditions equations ( Va ta co 2 phuong trinh dieu kien bien ) : ; 

. For determining the values at the interior mesh points we obtain the N-l equations ( De xac dinh cac gia 

tri cho cac diem trong , ta thu duoc N-l phuong trinh ) , then by replacing u '(x) and u "(x) ( Va thay the 

u'(x) va u"(x) ) : 



Ul := (k : t)-K 



1 < u k+l,t- u k-l,t> 



2 b 



U2:= (k : t)-»- 



u k + Lt- 2u k : t + u k-l,t 



We arrange this system of N+l equations in the form of matrix equation ( Sap xep he thong gom N+l 
phuong trinh nay ) . The matrix of it has N+l rows( Ma tran chinh co N+l hang ). The first row is fixed 
with the boundary condition at r = a ( Hang dau duoc xep cho dieu kien bien tai r = a ) . Obviously the last 
row is fixed with the boundary condition at r = b ( Hien nhien hang cuoi cung duoc xep cho dieu kien bien 
tai r = b ) . Now, we join these rows by listing them out , then construct the matrix symboled A . ( Lien 
ket cac hang nay lai , va xay dung nen ma tran A ) . 



The unknown values will be written as a vector ( cac gia tri chua biet se duoc viet dang vector ) 
u; : j=I..X 

and the right hand side of the equations is a column vector B ( va ve phai phuong trinh la 1 
vector cot B ) . Solving the matrix equation for u ( Giai phuong trinh ma tran tim nghiem u ) . Then we 
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f 100 



u(r : t) 



E 9 {r : t>:= 



.: 



+ 1 



find 



with 




and the expression of 



o^J,t) := 



C. NUMERICAL SOLUTION : 



Use Maple 9.5 

> restart : with(plots) : with(PDETools) : with(Linear Algebra) : 

m:=(100.+(t/To) A (l/10))*(t/To) A (2/5)/(100.+(t/To) A (l/2));To:=l; 

lame_cyl:=r A 2*diff(u(r,t),r$2)+r*diff(u(r,t),r)+m*u(r,t)=0; bound_con:=u(l,t)=l,u(2,t)=-l; a:=l; b:= 
N:=20; h:=(b-a)/N;R:=k->a+k*h; Ul:=(k,t)->(u[k+l,t]-u[k-l,t])/(2*h); U2:=(k,t)->(u[k+l,t]-2*u[k,t]- 
e[0]:=u[0,t]=rhs(bound_con[l]);e[N]:=u[N,t]=rhs(bound_con[2]); 

for k from 1 to N-l do e[k]:=eval(lame_cyl,{r=R(k),u(r,t)=u[k,t], diff(u(r,t),r)=Ul(k,t),diff(u(r,t),r$2^ 
row[0]:=[rhs(bound_con[l]),seq(0,j=l..N-l)]; 

row[l]:=[coeff(lhs(e[l]),u[0,t]),coeff(lhs(e[l]),u[l,t]),coeff(lhs(e[l]),u[2,t]),seq(0,j=l..N-3)]; 

for n from 2 to N-l do row[n]:=[seq(0,j=l..n-2),coeff(lhs(e[n]),u[n-l,t]),coeff(lhs(e[n]),u[n,t]),coeff(lh 

row[N]:=[seq(0j=l..N-l),rhs(bound_con[2])]; 

row_matrix:=[row[0],seq(row[n],n=l..N-2),row[N]];A:=(row_matrix); 

U:=Vector([seq(u[j,t],j=l..N)]);B:=Vector([rhs(bound_con[l]),seq(rhs(e[j])j=l..N-2),rhs(bound_coi 
print("Ham epsilon[theta](l,t)");plot3d(-U[N-l]/r,r=1..1.000001,t=0..100 );print("Ham epsilon[thetJ 
u(r,t):=U[N-l];;Ee:=0.5;;epsilon[theta](r,t):=u(r,t)/r;E[theta](r,t) := (100/(t/To) A .l+l)*Ee; sigma[th< 
sigma[theta](r,t)));sigma[b](t):=normal(subs(r=2,sigma[theta](r,t)));;;with(plottools):with(plots):plo 
style=[point,line],symbol=diamond,color=[red,black],thickness=[l],legend=[ sigma[bn](t) v , v sigma[b; 

Warning, the name changecoords has been redefined 
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m := 




To:= 1 




bound_con = n£l, t) = 1, u(2, t) = -1 



X:=:c 



h:=- 



R:=k-*a + kh 



l u k+l.t~ u k-l.t 
Ul:=(k : t)^- " 



U2 := (k ; t) ■ 



a k+l.t- Ju k.t + u k-l,t 



e 0 :=u 0,t =1 



e 20 : - u 20 ; t = - 1 



903 S61 
e l := ~ u 2.t- SS2u l,t + l" 1I 0,t + 



C2 := 495 u 3 t - 96S ip, t + 473 Uj t + - 



UOO. +t Jt uj t 
100. + Vt 

L GP 

V.100. + 1 Jt u 2 t 



1D0.+Vt 



"Ham epsilon[tlieta](l r t)" 
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"Ham epailon[tlieta](2 ; t)" 
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2.0000004 
2 0000006 r 
2.0000000 
100 2.000001 
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"HAMu(r : t)" 
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Live 




Disclaimer: While every effort has been made to validate the solutions in this worksheet, the authors are 
not responsible for any errors contained and are not liable for any damages resulting from the use of this 
material. 



Legal Notice : The copyright for this application is owned by the authors. Neither Maplesoft nor the 
authors are responsible for any errors contained within and are not liable for any damages resulting from 
the use of this material. This application is intended for non-commercial, non-profit use only. Contact the 
authors for permission if you wish to use this application in for-profit activities. 
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